%% =======================================================================
% TABLE 3 TOP
% =======================================================================
clear all; clc;
rng(11);
NN    = 1e6;

currentFile = mfilename( 'fullpath' );
[pathstr,~,~] = fileparts( currentFile );
addpath( fullfile( pathstr, '_tables' ) );

global cv_1stderror cv_2stderror

cv_1stderror = 1.00;
cv_2stderror = 1.96;

tt    = 1/15;
tx    = 0.01; 
tx1   = 0.06; 
teps  = 1.60; 
mu    = 0;

theta = mu*ones(NN,1)+sqrt(1/tt)*randn(NN,1);

v     = tx1/(tt+tx1);
vs    = tx/(tt+tx);
ty    = vs^2;

ei1   = sqrt(1/tx)*randn(NN,1);
xi1   = theta+ei1;

wx    = tx1/(tt+2*tx1+ty*teps);
wy    = ty*teps/(tt+2*tx1+ty*teps);

ei2   = sqrt(1/tx)*randn(NN,1);
eps   = sqrt(1/teps)*randn(NN,1);
xi2   = theta+ei2;
y     = (v/vs).*theta+(1/vs).*eps;

fi1  = v.*xi1+(1-v).*mu;
fi2  = wx.*xi1+wx.*xi2+wy.*y+(1-2*wx-wy).*mu;

error = theta-fi2;
mse   = mean(error.^2);
mse_p = (1-2*wx-wy*1)^2/tt+2*wx^2/tx1+wy^2/(vs^2*teps);   

cil  = fi2-cv_2stderror*sqrt(mse_p);
ciu  = fi2+cv_2stderror*sqrt(mse_p); 

count = zeros(NN,1);
for ii=1:NN
    outcome = theta(ii);
    if outcome>cil(ii) && outcome<ciu(ii)
        count(ii)=1;
    else
        count(ii)=0;
    end 
end
coverage_ratio_2std_model = sum(count)/length(count);

cil  = fi2-cv_1stderror*sqrt(mse_p);
ciu  = fi2+cv_1stderror*sqrt(mse_p); 

count = zeros(NN,1);
for ii=1:NN
    outcome = theta(ii);
    if outcome>cil(ii) && outcome<ciu(ii)
        count(ii)=1;
    else
        count(ii)=0;
    end 
end
coverage_ratio_1std_model = sum(count)/length(count);

confidence_bounds        = {'SPF-density-implied', 'Giordani-Soderlind-2003', 'Model-Implied'};
Percent95                = [0.85; 0.72 ;coverage_ratio_2std_model];
Percent66                = [0.59; 0.44 ;coverage_ratio_1std_model];
table_new                = table(confidence_bounds',Percent95, Percent66); 
writetable(table_new,'_tables/table_3_top.txt','Delimiter',' ')
table_new